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An EM-algorithm for dynamic SPECT tomography 

Heinz H. Bauschke 1 , Dominikus Noll 2 , Anna Celler 3 
and Jonathan M. Borwein 4 

1 Introduction 

Single Photon Emission Computed Tomography (SPECT) is a nuclear medicine diagno- 
sis technique which measures the three dimensional distribution of a radioactively labeled 
pharmaceutical injected in the body. As compared to standard imaging techniques like 
Computed Tomography (CT), the significance of SPECT lies in the fact that it reveals the 
function of the body rather than its structure. For example if a radio pharmaceutical is 
absorbed by an unhealthy tissue and rejected by healthy tissue, then SPECT will reveal the 
unhealthy tissue as a bright region. A related technique is Positron Emission Tomography 
(PET), [19]. 

A SPECT camera works by rotating around the patient an array of photo multipliers 
that detect gamma rays emitted by the patient. A collimator placed in front of the camera 
rejects rays that are not perpendicular to the camera face (see Figure 1). The images 
resulting in the camera are 2D projections of the original 3D activity distribution in the 
patient. Some devices use double or triple head cameras or even ring SPECT instruments 
to improve the number of detected counts and therefore the statistics. 

Current clinical applications of SPECT are based on the hypothesis that the injected 
radio activity in the organ of interest is stationary over the acquisition period (usually not 
more than 20 minutes). However, physiological processes in the body are usually dynamic, 
and some organs (kidney, heart) show a significant decay of activity due to wash-out. Being 
able to trace activity in space and time is therefore of importance and expected to signif- 
icantly enhance diagnostic possibilities. In particular, a dynamic SPECT reconstruction 
might be presented as a movie rather than a static image. 

With a significant decay of activity over the acquisition period, the filtered back pro- 
jection method (FBP) may no longer be reliably used, and new approaches have to be 
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developed. The purpose of the present paper is to present a probabilistic model for dy- 
namic SPECT whose numerical solution is based on an instance of the EM- algorithm. We 
discuss its numerical aspects along with those of other models as for instance presented in 
[13, 9]. 



2 Probabilistic Model 



In order to model dynamics in SPECT, we consider an instrument with a single camera 
head, available in most hospitals. The case of double or triple head cameras or even ring 
SPECT devices could be treated similarly. The camera rotates through 180°, with S stops, 
60 < S < 180, say, each stop taking 20/5 minutes. The camera plane is an array of photo 
multipliers of typical size 30cm x 40cm. Depending on its resolution, the camera is divided 
into M receptor bins, a typical choice being for instance some 6mm x 6mm, in which case 
we have M ~ 250 receptor elements. The camera rotates about a fixed axis, with a radius 
of 20cm - 30cm. 
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Figure 1. Photons radiating from the region of interest, a) misses the camera, 
b) is absorbed by the collimator, c) passes the collimator and hits the camera. 



The gamma rays originate from a region of interest of approximately 30 x 30 x 30cm 3 . This 
region of interest is divided into N little cubes called voxels. In our example, assuming a 
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spatial resolution of .5cm - 1cm leads to N = 30 3 up to N = 60 3 elements. During the 
kth stop the camera is at a fixed position with angle 6^ = (fc — 1) • ^y-, k = 1, ... ,5. 
What we are trying to reconstruct is the radio activity i(tk) of the isotope in voxels 

i = 1, . . . , N at times The projected data measured during each stop k present the 
activities in the receptor bins j = 1, . . . ,M. While the camera head is at position 9 k , only 
gammas traveling along a line with angle Ok are allowed to pass the collimator and hit some 
photo cell (see Figure 1). 

We assume that activity is constant during the time interval of a single stop, but is 
allowed to vary in time over the whole period of 20 minutes. Here activity Xik of voxel i 
during the fcth stop is proportional to the number of photons that leave the voxel during 
this time interval, radiating in any possible direction. 
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Figure 2. The coefficient c^k is the relative volume of voxel i lying within the 
beam Bjk connecting i and the receptor bin j during the camera position 0^. 



Let Yjk be the random variable which counts the number of events in the camera bin j 
during the stop fc. The physics of radio activity make it reasonable to assume that the Yjk 
are Poisson distributed. The data yjk collected in the camera bins at different positions 
form a sample for the Yjk- 

Let Cijk denote the probability that a photon leaves voxel i in direction Ok and hits the 
detector bin j while the camera is actually at position A:. The coefficients c^k are supposed 
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to be known. They are usually referred to as "the geometry" of the model, but besides 
geometry, they may as well include probabilistic factors like gamma ray scattering, and 
more importantly, attenuation (see [6, 20] for details on those). In a simplified model, 
where attenuation is ignored, Cij k is proportional to the volume of that part of voxel i which 
lies in the beam Bj k connecting i with the receptor bin j during the position k given by the 
angle 6 k , see Figure 2. Based on the geometry, the expected values of the Yj k are 

E(Y jk ) = ^2ci jk x ik >0. (1) 

For later use, let us fix some notation here. The coefficients C(j k give rise to two linear 
operators. We define C : R NS R MS by (Cx) jk = t c ijk x ik , and r : R NS — ► R NMS by 
(Tx)ij k = Cij k x ik . We observe that in practice T is usually injective, while C is typically 
not. 

3 Dynamic Tracing 

The critical part in modeling the dynamics of SPECT is the time behavior of activity 
Xi(t). If radioactive decay were the only significant part, a standard model of the form 
X((t) = Ai6~ Xt with known decay constant A > 0 would apply. Dynamic SPECT is however 
based on the hypothesis that due to a flow in the organ, each voxel i may have its own 
individual decay profile. Therefore, more sophisticated models for decay with rates varying 
in time and space are needed. Based on experimental evidence for fatty acid myocardial 
viability [10], the following parametric model of decay 

x t (t) = Aie'* + Bie-n* + G (2) 

with unknown parameters A { > 0,5; > 0,Ci > 0 and A* > 0, //< > 0 was used in [13, 3]. 
Discretizing (2) according to the 5 stops, t k = (k — 1)5/20, leads to a parametric model 
with 5 degrees of freedom for each z. Another approach suggested in the investigation of 
NMR relaxation data (cf. [21]) uses a model of the form 

Xi (t)= f a;(A)e- A '<ZA. (3) 
Jo 

Discretizing the time axis at the S stops t k now leads to a sum of exponentials with S 
degrees of freedom. 

A third approach which we propose here is to not insist on any specific form of the decay 
curves. That is, in each position i we allow for an activity profile with 5 degrees of freedom 
x iky k = 1, . . . , 5, assuming only that activity decays in time: 

> %i2 ^ £i3 > ■ ■ • > Xis and Xi$ > 0. (4) 

As compared to (2), this increases the flexibility of modeling decay, but increasing the de- 
grees of freedom bares the risk of producing an underfitted model. Fortunately, experiments 
indicate (see section 8) that the model works well, and that its obvious algorithmic and nu- 
meric advantages over (2) should be exploited. To introduce some standing notation, we 
denote the set of x = {x ik } satisfying (4) by Q,. 
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4 Maximum Likelihood 



Let us now derive the maximum likelihood estimation for the unknown parameters x € Q 
based on the data y. Let g{y\ x) be the probability density of the measurements Y given 
the activity curves x. Then we consider 

(ml ) maximize £( lo g£(y; x )) 

^ x ' subject to x e £1 

i.e., we pick those parameters x € Q which were the ones most likely to have produced 
the data y — {yjk}- As the Yjk are Poisson variables with means X^ili c ijk%ik-> up to some 
constant the negative log-likelihood function — \ogg(y\x) equals 

MSN N 

<P( X ) — Yl 2 ( 2 Ci i kXik ~ y i k log ( S Ci J kXik ) ) * 

j=\ k=l i=l t=l 

The problem of minimizing <p(x) subject to the constraints (4) is feasible and admits optimal 
solutions, since the objective clearly tends to infinity if at least one of the x ik — > oo on the 
feasible set fL On the other hand, solutions are typically not unique. Any solution x e ft 
must satisfy the following Kuhn- Tucker conditions (cf. [7]): There exist multipliers A^ > 0, 
i = 1, . . . , N and k = 1, . . . , S such that for each i — 1, . . . , N, 

N 



j=l x'=\ 

M N 



(/vTO 



j=l t'=l 

M N 



Tis-i ~ X ( C U"S-12/?S~1 / X ^'i^-l^'^-l) +Ks-2 -A<5_i = 0 

j=l i'=l 

j=l i'=l 

£ik > Aijt^jt - a;tjb+i) = 0 fc = l,... ) 5-l, A^s = 0 

where we used the abbreviation 

M 

r ik := Yl Ci i k > °* 
i=i 

Clearly any solution x has £) £ cy*^* > 0 for fixed j, fc. Notice also that any two solutions 
x,i E fi must satisfy Yli c ijh^ik = ^(CijkXik for every fixed j, fc, since the functions 77 — * 
77 — yjk log?; are strictly convex for yjk > 0. As the operator given by (Cx)jk = ]Cjli c ijk x ik 
is typically not injective (MLi) may have multiple solutions. (Notice that ensuring that C 
is only moderately defective is of relevance for the numerics and may be influenced by the 
discretization we choose; cf. [3]). 
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5 The EM-algorithm 



The EM-algorithm is an iterative procedure to calculate maximum likelihood estimates. It 
first made its appearance in 1976 in a paper by Dempster, Laird and Rubin [8], and its 
original intention was maximum likelihood estimation with incomplete data. Since then, 
the EM-algorithm has been applied in a much wider context, including situations in which 
the incomplete and the complete data space are defined in a somewhat artificial way. 

In our present situation, the maximum likelihood problem (ML\) involving the joint 
probability distribution g(y\ x) of Y represents the incomplete data space. To introduce a 
complete data space we consider the random variables Zijk which present that part of the 
activity in voxel i that radiates towards the receptor bin j during stop k. These are Poisson 
random variables with mean c^rr^ and joint probability density f{z\ x) depending on the 
parameter x E ft. 

The EM-algorithm is now the following. By induction define a sequence of parameter 
estimates x a = {xf h } 6 ft, a = 1, 2, . . . according to the following rules: 

1. E-step: Given x a = {xf k } € ft, and the data y, calculate the conditional expectation 
z a = E(Z\Y = y,x a ). 

2. M-step: Calculate the new estimate x a+l € ft by maximizing E(logf(z a ]x) \ y\x a ) 
over x £ ft. 

The E-step is fully explicit here. Namely, given the data y and the current x a , we find that 

z ijk ~ Vjh • W 



^ C i'jkX?'k 



i'=\ 

This may in fact be obtained from the following: 

Lemma 1. Let X\, . . . . X n be independent Poisson distributed random variables with mean 
'E(Xk) = /4t- Then the conditional expectation 

E((X X , . . . , X n ) | X x + . . . + X n = d) =: (xj , . . . , x n ) 

is obtained as 

, Vk 

Xk = a • 



□ 

For the M-step, we maximize the log-likelihood function logf(z a+1 \x) over the x € Q. 
Up to constant terms, that entails minimizing 



N M 



#0 -= ^2 X £ (^ifc^fc - Z ijk log(^'^ifc)) 



t=l j=\ k=\ 
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over all x € fi. Obviously, this function is of the form ift = £V ifii with each 



M S 



depending only on the variable = (xn, . . . , #is) € fit, where 17^ are the x l satisfying (4) 
for fixed i. Notice that the data y and the previous iterate x a do not enter directly into the 
calculation of x a+1 . Rather it is z a which makes connection with the previous step. 
We have 



--l^[ C M-—)- T *-—' 

ik . * x ^ik ' ^ik 



with Tjfc as in (5) and 



M 

j=i 

We now consider the maximum likelihood problem for the complete data space, which 
perforce splits into N problems of size 5: 

(ML ) maximize E(logf(z a : x) \y;x a ) 
^ 2) subject to x e ft 

Setting ipik{x l ) := —Xik + x^+i, and with ibi defined as before, the ith problem becomes: 

(MT , minimize ?pi(x l ) 

{ < ML2 > i > subject to ^ifc(^) < 0 (fc = L . . . , S - 1) 

Notice that in contrast with the original maximum likelihood problem (MLi), we do not 
control the constraint Xis > 0 here, since it is built in through the objective due to a? s > 0. 
In fact, for fixed z, k, consider j such that both c^k > 0 and zf jh > 0. Then the corresponding 
term in ipi(x l ) works as a barrier function which forces the variable Xik to take strictly 
positive values. For the same reason, it is now reasonably clear that the optimal solutions 
:r l,a+1 for (ML2,;) a re unique, since the ipi are strictly convex and coercive. 

The Kuhn-Tucker conditions (cf. [7]) for the optimal solutions x t,a+} of (ML2,i) imply 
the existence of multipliers A^. > 0, k = 1, . . . , S — 1 such that: 

n ~\ a - fi 



r i2 -^r +AS -AS = 0 



(KT2,i) 



X i,S-l 



S 



x^ 1 > xfi^i Kk( x ik l ~~ x Tk+i) — 0 /c — 1, 2, . . . , S — 1. 
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Based on the optimality conditions (KT 2) i), we may now enter a detailed analysis of the 
present EM algorithm. 

6 Alternating Projections 

We show that the EM-algorithm presented in the previous Section may be interpreted as an 
alternating projection procedure (in the sense of von Neumann) provided that the Euclidean 
metric is replaced with the Kullback-Leibler distance. 

Let n = NS and m = MS, p = NMS. The Kullback-Leibler distance of a, b € is 
defined as 

v 

d(a, b) = ai lo S( a */ b i) ~ a i + 

Now d has properties which resemble those of a metric, but obviously lacks symmetry and 
is only defined for a > 0 and b > 0. In analogy with the orthogonal projection we define 
projection operators and PJJ~ associated, with d. Given a closed convex subset A of 
and a point x G R+, the forward resp. backward projections of x onto A are defined as 

and 

PZ(z) = argmin ae ^d(x,a). 

With these notions, the M-step presented in Section 5 turns out to be a Kullback-Leibler 
backward projection. 

Proposition 1. Let A := {v € R NMS : = c ijk x ih for some x e fi}. Given 2° e R NMS 7 
let t> a+1 6e £/ie backward projection of z a onto A, that is 

^a+l = p-(^). 

Tften tfie solution x a + l of the M-step satisfies vf^ 1 = (H^xf^ 1 , for short, i> a+1 = Fx a+l . 

Proof. Indeed, for the proof we simply observe that with Vijk = Qjjt^ii, for short, v = Tx, 
the objective function satisfies tp(x) = ]Cijfc( v tjJb — zfjk l°E v ijk)i which up to a constant term 
equals d(z a ,v). Minimizing ip(x) over x G H is therefore equivalent to minimizing v) 
over v £ A. Existence and unicity of the projection is guaranteed here since z a > 0 and the 
set A contains points in R++ s . Notice also that A is closed as it is the linear image of the 
closed convex and polyhedral set fi under T, and is therefore itself polyhedral. □ 
Based on this observation, we now ask whether similarly the E-step may be interpreted 
as a projection onto a convex set. 

Proposition 2. Let B :={z € R NMS : z > 0, z ijk = y jk Vj, k}. £e* x Q be the current 

EM-iterate, and let v a = rV*. JTien the conditional expectation z a defined by the E-step 
(6) is the Kullback-Leibler projection of v a onto B in the forward and the backward sense: 

z° = P^(^) = Pf(v a )- 
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Proof. First consider the case of the backward projection. Clearly z = Pjf(v a ) exists and 
satisfies the Kuhn-Tucker conditions: There exist multipliers Xj^ such that 

-^ + 1 + A ijb = 0 Vi,j,fc 

Zijk 

in tandem with the constraints 

N 

t=l 

Summing over i for fixed j, k gives (1 + Xjk)yjk = Si reac lily implies formula (6). 

Next consider the case of the forward projection z = Pj3*(v a ). Here the Kuhn-Tucker 
conditions provide multipliers Aj*. satisfying 

l°g(zijk/v? jk ) + \jk =0 Vi, j, k 

and with constraints as above. Taking exponentials and summing over i for fixed j, k gives 
the same type of relation and again leads to formula (6). □ 
It remains to observe that any limit point of the sequence x a must be a Kuhn-Tucker 
point for the problem (ML\). This is a consequence of a general fact (see [22]), but could 
equally well be checked using the conditions (KT\) in tandem with (6), (7) and the (A^i). 
Indeed, convergence of the iterates v a } z a implies convergence of x a , <j a . Formula (KKTi) 
implies 

k 

^k = Y,(T ij -a« j /x° j l 
j=i 

so the multipliers converge as well. Passing to the limit a — ► oo therefore implies the 
corresponding conditions (AT\), and it remains to check that part of the Kuhn-Tucker 
conditions (AT\) which concerns the constraints X{g > 0. But clearly, if any such constraint 
is active, Xi$ — 0, the corresponding conditions in (KT\) are satisfied. In case Xi$ > 0, the 
multiplier Xis = 0 will do. 

In consequence, we may now state the principal convergence result for the present Pois- 
son EM-algorithm. 

Theorem 1. The sequences v a = Tx a and z Q generated by the Poisson EM-algorithm based 
on the constraint set (4) converge to limit points v a — > v € A, z a — » z € B. Here v is of 
the form v = Fx, with x € ft a solution to the maximum likelihood problem (ML\). 

Proof. First observe that the problem (ML\) is convex and, due to the structure of the 
constraint set ft, admits optimal solutions. We claim that the sequences of alternating 
Kullback-Leibler projections 

z a = Pb(v°) = Bg(v a ) v a+1 = PZ(z a ) 
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converge to certain limit points v a — > v G A, z a — ► z € 5. In the case where ^4fl5 ^ 0, the 
alternating sequence converges to a common limit point v — z £ AC\ B. This is a known 
fact proved by Bregman (see e.g.[l, Thm. 8.1]). In case A fl B = 0, which is the more 
realistic one, the convergence proof is more complicated, and only special cases have been 
considered. See for instance [11, 12, 19, 4], where the case of positivity constraints x ik > 0 
was discussed. Our present case involving the constraint set (4) may be settled using the 
same type of reasoning. We skip over the details, which are tedious. 

Finally, notice that the point x satisfying v = Tx is uniquely determined by the injec- 
tivity of T, and by closedness lies in Q. Inspecting the Kuhn-Tucker conditions for (MLi), 
we see that x must be a Kuhn-Tucker point, and by convexity therefore solves (ML\). □ 

Remark. Consequently, the EM-algorithm is best understood in the space of ^-variables, 
where it may be recast as an alternating Kullback-Leibler projection onto closed convex 
sets A and B. In fact, B is an affine subspace, restricted to while A = T(fi) depends on 
the parameter space 0, and may therefore be quite complicated. The whole procedure may 
indeed be generalized to any comparable situation with a closed convex parameter space fL 

Remark. Replacing the constraint set Q by fl representing the nonlinear conditions (2), we 
obtain a non-convex set A = T(Q). The EM-algorithm could then still be formulated as an 
alternating projection between A and 5, but with the obvious problems when projecting 
onto non-convex sets. 

Remark. The case considered by Iusem [11] and others becomes the case where we impose 
only positivity conditions, Cl = {x : rr^ > 0}, A = T(Q). In this case the backward 
projection onto A is particularly pleasant to calculate and leads to an explicit formula. 



7 The Least Squares Approach 

In [9], the authors propose a different approach to the dynamic SPECT problem, which is 
based on the dynamic model (2) and is solved using nonlinear least squares: 

minimize \\Cx — y|| 
(NLS) subject to x ik = Aie~ Xiik + Bie"^ itk + d 

Ai > 0,B { > 0,Ci > 0,A; > 0,^ > 0 

Here the operator C : R n — ► lR m is defined as before: 

N 

Problem (NLS) was found difficult to treat in practice, due to the highly nonlinear depen- 
dence of the x^ on the parameters (cf. [13, 3]). Following the ideas presented in Section 
3, we replace the parametric model (2) by (4). This leads to a linear least squares problem 
with inequality constraints 

(LLS) minimize ~ y ll 

^ ' subject to x £ f) 
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where Q is as before the set of x satisfying (4). Our aim now is to present a probabilistic 
model for both (NLS) and (LLS). This will also provide clues to possible numerical 
approaches based on a version of the EM-algorithm. 

Following the pathway in Section 5, we let the random variables Z^ and Yjk be defined 
accordingly, but now assume that the Z^*. are independent and normal variables with 
unknown mean CijkXik and common variance a 2 > 0. Consequently, the Yjk are independent 
normal variables with mean 5^ i==1 Cijk x ik and variance No 2 . 

With these agreements, (LLS) is just the maximum-likelihood estimation problem 



where g(y\ x) denotes the density of the normal law Y with mean Cx and covariance matrix 
Na 2 I m . Similarly, problem (NLS) is the maximum likelihood problem with the constraint 
set Q, given by (2). 

Returning to the outline of Section 5, we apply the EM-algorithm to the present situ- 
ation, letting Z resp. Y represent the complete resp. incomplete data spaces. As in the 
Poisson case, this generates sequences x a and z a , a = 1, 2, . . . according to the following 
rules: 

1. E-step. Given the current iterate x a € fi, calculate the conditional expectation 
z a = E(Z\Y = y;x a ). 

2. M-step. Calculate the new iterate x a+1 € SI by maximizing E(logf(z a ]x)\y]x a ) 
over x e fi. 

Here f(z\x) denotes the normal law with mean Fx and covariance matrix o 2 I v , which is 
the joint density of Z. 

Proposition 3. Let A = {v € R p : v = for some x E ft}. Then, given the result z a of 
the previous E-step, the next M-step reduces to calculating the orthogonal projection v a+l of 
z a onto A, and taking x a + l to satisfy v a + l = rx Q+1 . 

Proof. Indeed, the negative log-likelihood function used for the M-step up to constant 
terms equals \\z — rx\\ 2 /(Na 2 ), and is to be minimized over x E ft. This establishes the 
statement. □ 
Let us now pass to the E-step, which we would like to reveal as an orthogonal projection. 
This involves calculating the conditional expectation E(Z\ %ijk — Ujk]z a )- Due to 
independence of Z, this may be obtained from the following: 

Lemma 2. Let X = (X\, . . . y X n ) be a vector of independent normal variables with mean 
E(X{) — fa and variance a 2 = 1. Then the conditional expectation x = E(X \ X\ + . . . + 
X n = y) satisfies 



More precisely } x is the orthogonal projection of . ,/i n ) onto the set of (x\, . . . ,x n ) 

satisfying x\ + . . . + x n = y. 



{ML,) 



maximize E( \ogg(y; x)) 
subject to x € 
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Proof. In consequence of [18, Thm. 2.1(viii)], a normal vector 



gives the conditional expectation 

E(b | a) ~ M((3 + W*U-\a - a), V - W*U~ l W). 

We apply this in the case b = (X\, . . . 5 X n ), a = X\ + . . . + X n , which gives the stated 
formula. □ 
With this observation we find that the formula replacing (6) in the case of the normal 
variables is 

1 1 N 

z tjk = x yjk + CijkX * k ~ N S Ci 'i kX "k', (8) 

t'=i 

and the E-step is in fact realized as an orthogonal projection: 

Proposition 4. Let B = {z £ R p : YliLi z ijk — Vjh Vj, k). Let x a be the current iterate 
generated by the previous M-step, and let v a = Fx a . Then the result z a of the next E-step 
is the orthogonal projection of v a onto B. □ 
We are now in the position to state the convergence result for the EM-algorithm in the 
case of normal variables. It is based on von Neumann's classical method of alternating 
projections. 

Theorem 2. Let A^B be defined as above, and let z a , v a = Fx a be the sequences generated 
by the Gaussian EM-algorithm. Then z a , v a converge to certain limit points v a — ► v £ 
A, z a — ► z € B. Here v is of the form form v = Fx, with x € fi a solution to the least 
squares problem (LLS). 

Proof. In the case where A fl B ^ 0, it is well-known that the sequence of alternating 
projections converges to a common limit point v — z (E A Pi B, which by definition of A is 
then of the form Fx. 

The more involved case occurs when A Pi B ~ 0. Following [2], a dichotomy appears: 
Either the alternating sequences converge to limit points v a ^ v € A, z Q — > z € B : with v, z 
realizing the distance between A, B, or there are no points realizing this distance, in which 
case the sequences tend to infinity, \\v a \\ — ► oo, \\z a \\ — ► oo, with \\v a - z a \\ approaching the 
distance between A y B. We will argue that in our situation the second case is impossible. 

Indeed, the distance d(A,B) not being attained implies the existence of a common 
asymptotic direction for A, B: there exist iu ^ 0 and a G A, b € B having a + R+w C A, 
b + R+w C B. The second inclusion implies Y^i w ijk = 0 for fixed j, k, while the first gives 
w = Fx for some x € fi. The c ijk being nonriegative, this is only possible when x = 0, and 
hence w = 0. 

We deduce that v a — ► v = Fx and z a — ► z. It is now routine to check that x is a 
Kuhn- Tucker point for the original problem (LLS). □ 
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Remark. Replacing Q by 0 = {x : x^ > 0} again gives an important special case. The 
orthogonal projection onto A = T(Cl) may be calculated explicitly, and in tandem with 
(8) leads to an explicit formula for updating x a . We then have an iterative method for 
calculating nonnegative least squares solutions. One may compare this to other iterative 
methods for calculating nonnegative least squares solutions like ISRA (see [16]). 

Let us analyze the M-step a little further. Recall that up to a constant additive term 
and a constant factor, the negative log-likelihood function — log f(z a .x) equals 



2 - ' l 

and is to be minimized over Q. Setting 

M M 

nk == 4* > 0 and 6ik := > °> ( 9 ) 

we have 

ij){x) = - {lik%ik ~ %bik%ik) + constant terms. 

The Kuhn- Tucker conditions therefore imply the existence of Lagrange multipliers A^ > 0, 
i = 1, . . . , N and k = 1, . . . , S satisfying 

Ii2%i2 -fii2 + - A i2 =0 

(LKTi) 

JiS-lXiS-l -$is-\ +Ks-2 -A,-5_i = 0 

lis%is -$is +Afs-i -\is = 0 

Xik>*ik+u ^ik&ik - ^tJt+i) = 0. for k = 1,... 5-1, ^5 > 0, A^s = 0. 

So analogously to the Poisson case, the original problem of size NS splits into N problems 
of size 5. Notice that in contrast with the Poisson case (ML\^), we do have to control the 
constraints Xi$ > 0 here, which leads to the additional multipliers A^. However, as we shall 
see, these constraints will generally be inactive, and the multipliers will vanish. 



8 Our Numerical Approach 

In this section we shall discuss the practical aspects of both EM- algorithms. Clearly the 
algorithms may be expected to be slow, since alternating projections are known to converge 
no better than with a linear rate. Nevertheless, the additional numerical stability gained 
may often justify using an EM-scheme, in particular, if some speed is recovered e.g. by 
parallelizing the M-step. The crucial question to be addressed before proposing this scheme 
is whether the M-step, which is intrinsically more complicated than for the stationary case 
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([11, 4]), allows for a fast numerical solution. Let us therefore analyze the M-steps of both 
the Poisson and the Gaussian models. 

The first surprising observation is that both M-steps, although coming from completely 
different out-sets, lead to exactly the same procedure. Recall the Kuhn-Tucker conditions 
(7\T 2? i) for the Poisson M-step presented in Section 5. To simplify notation we omit the 
indices i and a which are fixed, writing r* := and Ok := and similarly x k := xf h for 
the solution and the multipliers. We then have the following: 

Lemma 3. For fixed i the problem (ML2,i) has a unique solution (x\, . . . . #$) satisfying 
(4). This solution is of the following form 

X\ = . . . = x ri > # ri _|_i = • . . = x T2 > . . . > x rt _ 1 +i — . . . = x rt (10) 

for appropriate 0 = ro < r\ < . . . < r t = S, with 



„ _ _ (j) _ frrj.i+l + ■■■ +Vrj 



(11) 



*/ie being strictly decreasing. In particular x\ = x^ > ^ and x$ — x^ < 

Proof. Uniqueness of the solution of (ML 2 ,i) follows from strict convexity of the objective. 
Clearly any solution, since it satisfies (4), has the form (10), and it only remains to establish 

(ii). 

The complementarity condition in (AT^i) implies that A r> = 0 for all j = 1, . . . 
while the remaining multipliers may be strictly positive. Now summing the Kuhn-Tucker 
equations (/\T 2jt ) in each block rj_i + 1, . . . , rj separately gives formula (11). □ 

Remark. This result is interesting since it tells us that given any sequences a k > 0, r fc > 0, 
k = 1, . . . , 5, there exists a unique subdivision (10) such that the sequence x^\ j = 1, . . . , t 
defined through (11) is strictly decreasing. (Taking for instance r k = 1 we obtain the 
amusing observation that, given any sequence > 0, there exist a unique subdivision (10) 
such that the arithmetic means of the Oi over each block are strictly decreasing and satisfy 
the boundary conditions x^ > o\ and x^ < 

Let us now pass to the Gaussian M-step obtained in Section 7. Here we have the same 
observation, which, again on suppressing the indices i,a, is the following: 

Lemma 4. The solution (x\ ..... x$) of the Kuhn-Tucker conditions (LKT{) is unique and 
of the form (10), possibly with a different t. It admits a representation of the form (11), 
with 6i replacing and 7,- replacing r,-. The sequence x^ is strictly decreasing, and again 
the limit conditions X\ = x^ > x$ = < ~ are satisfied. 

Proof. Indeed, starting out with the blocks (10), possibly with another t, we find again 
that the multipliers \ Tj must vanish. Again summing the Kuhn-Tucker conditions (LKT{) 
in each block would suffice, provided we knew that the final multiplier A 5 belonging to the 
constraint x$ > 0 vanished. Clearly this is the case if x$ > 0, which is consequently what 
remains to be checked. 



14 





broken line - model (2) 
continuous lino - modal (4) 






I 



Comparison of model* (2) and (4). with 60 stops 



10 20 30 40 SO 60 70 80 

Comparison of model* (2) and (4). with 80 stop* 



broken lina - modal (2) 
continuous Una - modal (4) 
dotted line - adapting (2) to (4) 



20 30 



50 60 



Comparison of (2). (4) and smoothing (4) via <2) 
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Comparison of (2). (4) and smoothing (4) via (2) 



Figure 3. Comparisons of models (2) and (4) using randomly generated data. In the lower line an 
a posteriori fit of (2) to the data of (4) shows little difference with model (2). 



Assume to the contrary that the last block in (10) is zero: = £ r ,_i+i = . . . = %s = 0. 
Then the Kuhn- Tucker conditions give 

5 

which in view of As > 0 and 6k > 0 is only possible when <5 r< _ 1+ i = . . . = 6s ~ 0. However, 
by construction (9), we have 6i > 0, so \$ — 0, x$ > 0, and the result follows. □ 
The observation that both M-steps are essentially equivalent, although with different 
data crjt > 0, Tjt > 0 versus 6k > 0, 7^ > 0 gives us choice on how to perform the M-step. In 
fact, the Gaussian M-step uses nonnegative least squares, and for a moderate size S works 
faster than the Poisson M-step (ML2,i). On the other hand, for a really large 5, if required, 
the Poisson M-step could more conveniently be solved by an interior point method. In fact, 
a logarithmic barrier term for the constraints (4) leads to the objective: 

S S-l 

One might perform only a few steps towards minimizing cp^ for a fixed fi > 0, before passing 
to another E-step. Increasing the number of steps and reducing /i should then be controlled 



15 



by testing decrease of the likelihood function in {ML\) resp. in (LLS), which for theoretical 
reasons (see [8]) is known to decrease at each iteration if an exact EM-step is performed. 

Experiment. Figure 3 shows typical results of the Gaussian M-step for randomly gener- 
ated data on the basis of 60 (left column) resp. 80 (right column) camera positions. The 
plots in the top line show comparisons of typical decay curves produced by the models (2) 
(smooth curves) versus (4) (steplike curves). 

The lower pictures indicate that typically no information is lost on replacing (2) with 
(4). Namely, the smooth dotted lines show model (2) curves which a posteriori have been 
fitted to the steplike functions (4). Typically, these show little difference from the direct 
fitting of (2) (smooth broken curves as before), so (2), if desired, could be retrieved from 
(4). As (4) has obvious numerical advantages, the EM-algorithm should in fact be built on 
(4). Fitting a model (2) decay curve to the step functions (4) may be deferred to the end 
of the procedure. 

9 Cyclic Projections 

The fact that the Kullback-Leibler resp. orthogonal projections onto the set A cannot be 
calculated explicitly is clearly a drawback of our present approach. This may be overcome 
by splitting the alternating projections based on the EM-algorithm into three or even more 
projection steps which are easier and may be performed in a cyclic ordering. To be more 
precise let us consider the sets 

A\ {v 6 R p : v = Tx, Xn > x i2 , X& > a;^, . . .} 

and 

A 2 := {v G R p : v = Fx, x i2 > x i3 , x iA > . . .}. 

Then A — A\ H A 2 , and instead of projecting onto A, we wish to use the projections onto 
the Ai, which as we shall see are easier to perform. 

Inspecting a simple case where the sets A\, A 2 , B are the edges of an equilateral triangle 
show that projecting cyclically onto the sets is not exactly what we want. The are, however, 
various possible ways in which the projected information onto the three sets may be used to 
approach the points a G A, b 6 B realizing the distance between the two sets. We propose 
the following scheme, which applies to both the Poisson and the Gaussian case with the 
corresponding interpretations: 

1. Given an iterate x a resp. v a = Tx a , do an E-step using (6) resp. (8). The result is 

2. Replace the M-step by the following: Project z a onto A{, which gives v l = P^iz**), 
i = 1,2. To form the new iterate x a+1 take v a + l = (v x + v 2 )/2, v a + l = r(x a+1 ). 

3. At the end of the procedure do a few correction steps by projecting v a onto the true 
set A. 
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Let us see why this scheme, leading to explicit formulas, is expected to be faster than 
the original EM-algorithm. The crucial observation is the following 

Lemma 5. With the notation (7) and (5): Let v l = be the Kullback-Leibler back- 

ward projection of z a onto A\ . Then v l = Fx* and x l is given by the following alternative: 
Either a iy 2k-\/Ti,2k-i < <^,2fc/ r i\2Jfe; in which case 

or <Ti£k-ifn,2k-i > 0i,2k/Ti,2k, in which case we have 

X \,2k-\ - ^t,2*r-l/Tt,2Jfe-l, an d ^i,2Jfc = Vi,2k/Ti,2k- 

A similar formula is obtained for the backward projection v 2 = f%(^ a )j ^ 2 = T^ 2 - □ 

10 Conclusion 

We have presented two versions of an EM algorithm for dynamic SPECT based on a Poisson 
resp. a normal distribution. The EM algorithm is recognized as an alternating projection 
scheme either in the sense of von Neumann, or with respect to the Kullback-Leibler distance. 
Our simulations indicate that a parametric model based on experiments in myocardial via- 
bility studies (2) used in previous experiments [9] should be replaced by a less biased model 
(4), which in addition has computational advantages. Modified versions of our approach 
(in the spirit of [17]) accounting for measurement noise in real clinical data could easily be 
formulated. Related regularization techniques are discussed in [14, 15]. 
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Abstract — In this paper we present two variants of the EM 
algorithm for dynamic SPECT imaging. A version based on com- 
partmental modeling which fits a sum of exponentials and a more 
general approach allowing for arbitrary decaying activities. The 
underlying probabilistic models are discussed and the incomplete 
and complete data spaces are shown to be physically meaningful. 

We indicate that the second method, leading to a convex pro- 
gram in the M step, is easier to treat numerically and we present 
a possible numerical approach. Some preliminary numerical tests 
indicating the feasibility of the method are included. 



I. Introduction 

SINGLE photon emission computed tomography (SPECT) 
is a nuclear medicine diagnosis technique which measures 
the three dimensional (3-D) distribution of a radioactively 
labeled pharmaceutical injected in the body. As compared to 
standard imaging techniques such as computed tomography 
(CT), the significance of SPECT lies in the fact that it reveals 
the function of the body rather than its structure. For example, 
if a radio pharmaceutical is absorbed by an unhealthy tissue 
and rejected by healthy tissue, then SPECT will reveal the 
unhealthy tissue as a bright region. A related technique is 
positron emission tomography (PET), [34]. 

A SPECT camera works by rotating around the patient a 
scintillation detector (camera head) that records gamma rays 
emitted by the patient. A collimator placed in front of the 
detector rejects rays that are not perpendicular to the camera 
face (see Fig. 1). The images resulting in the camera are 
two-dimensional 2-D projections of the original 3-D activity 
distribution in the patient. Some devices use double or triple 
head cameras, or even ring SPECT instruments, to improve 
the number of detected counts and therefore the statistics. 

Current clinical applications of SPECT are based on the 
hypothesis that the injected radio activity in the organ of 
interest is stationary over the acquisition period (usually about 
20 min). However, physiological processes in the body are 
usually dynamic and some organs (kidney and heart) show 
a significant decay of activity due to washout. Being able to 
trace activity in space and time is therefore of importance and 
is expected to significantly enhance diagnostic possibilities. A 
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Fig. 1 . Photons radiating from the region of interest, (a) Misses the camera, 
(b) Absorbed by the collimator, (c) Passes the collimator and hits the camera. 



dynamic SPECT reconstruction, if numerically and algorith- 
mically possible, might be presented as a movie rather than 
a static image. 

With a significant decay of tracer activity over the ac- 
quisition period, the filtered back projection method (FBP) 
may no longer be reliably used, and new approaches must be 
developed. The purpose of the present paper is to introduce 
a probabilistic model for dynamic SPECT whose numerical 
solution is based on an instance of the EM algorithm. We 
discuss its numerical aspects along with those of other models 
as, for instance, presented in [20] 

Let us mention that dynamic SPECT is currently already 
possible with a ring SPECT instrument, rarely available in 
hospitals, or with a triple head camera doing fast camera rota- 
tions. Here, the time-varying activity is handled by acquiring 
a number of fast sequential SPECT studies and fitting curves 
on a voxel-by-voxel basis to the resulting SPECT images (cf. 
[30]). This approach suffers from poor data statistics and is not 
compatible with the majority of the transmission scan-based 
attenuation-correction methods (cf. [32], [31], [9], [18], [19]). 
Our investigation aims at a dynamic method which works even 
for single-head cameras doing slow rotations adapted to the 
time scale of the tracer dynamics. 

Clinical applications of dynamic SPECT which could be 
realized in the near future include renal studies using the 
dynamic tracers Tc-99mMAG3 and Tc-99mDTPA used in 
planar renal imaging, myocardial viability studies using 1-123 
labeled fatty acids, or brain blood-flow imaging using Xe-133, 
Xel27, or 1-123, (cf. [15], [23]). The biological half lives of 
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any of these dynamic tracers is within the time scale of a 
typical scan of 20-25 min duration. 

II. Probabilistic model 

In order to model dynamics in SPECT we consider an 
instrument with a single camera head, available in most 
hospitals. The case of double- or triple- head cameras or even- 
ring SPECT devices could be treated similarly. The camera 
rotates through 180° with 5 stops, 60 < S < 120, with 
each stop taking T/S minutes. T is the total acquisition 
time, usually some 20 min. The camera plane is of typical 
size 30 x 40 cm. Depending on its resolution, the camera is 
divided into M receptor bins (pixels), a typical choice being 
some 6 x 6 mm, in which case we have M « 3000 receptor 
elements. The camera rotates about a fixed axis with a radius 
of 20-30 cm. 

The gamma rays originate from a region of interest of 
approximately 30 x 30 x 30 cm 3 . This region of interest 
is divided into N little cubes called voxels. In our example, 
assuming a spatial resolution of 0.5-1 cm leads to N = 30 3 
up to iV = 60 3 elements. During the fcth stop the camera 
is at a fixed position with angle 9k = {k - 1) • 180°/5 J 
k = 1, ... ,5. What we are trying to reconstruct is the radio 
activity Xik := Xi(tk) of the isotope in voxels i = 1, . . . ,N 
at times i*. The projected data yjk measured during each stop 
k present the activities in the receptor bins j = 1,...,M. 
Ideally, while the camera head is at position k only gammas 
traveling along a line with angle are allowed to pass the 
collimator and hit some photo cell (see Fig. 1). 

We assume that activity is constant during the time interval 
of a single stop, but is allowed to vary in time over the whole 
acquisition period. Here, activity x^k of voxel % during the 
fcth stop is proportional to the number of photons that leave 
the voxel during this time interval, radiating in any possible 
direction. 

Let Yjk be the random variable which counts the number 
of events in the camera bin j during the stop k. The physics 
of radio activity make it reasonable to assume that the Yjk are 
Poisson distributed. The data yjk collected in the camera bins 
at different positions form a sample for the Yjk- 

Let Cijk denote the probability that a photon leaves voxel 
i in direction 9k and hits the detector bin j while the camera 
is actually at position k. The coefficients Cijk are supposed 
to be known. They are usually referred to as the geometry 
of the model but, in addition to geometry, they may include 
probabilistic effects such as gamma ray scattering and, more 
importantly, attenuation and collimator blurring (see [1], [9], 
[36], [33]). 

While including the collimator response function into the 
model parameters Cijk is routinely possible, correcting for 
attenuation causes a major problem which, in practice, is 
handled in two ways. The first is to calculate the attenuation 
map beforehand, either by doing a transmission scan in parallel 
with the SPECT acquisition (see for instance [9]) or by 
estimating attenuation mathematically. For the latter approach 
see the method proposed in [24] and reported to be practical 
in [35]. A second way, at least mathematically possible, 
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Fig. 2. The coefficient c^is the relative volume of voxel i lying within the 
beam Bj k connecting i and the receptor bin j during the camera position 9 k . 

consists in estimating the unknown attenuation and source 
simultaneously as proposed, for instance, in [6] and, more 
recently, in [13]. 

The way attenuation is estimated does not have any influ- 
ence on our dynamic protocols. The coefficients Cijk may, 
in principle, even comprise scattering, although the number 
of nonzero Cijk then increases dramatically, leading to a 
numerically difficult problem. 

In a simplified model where collimator blurring, attenuation, 
and scatter are ignored Cijk is proportional to the volume of 
that part of voxel i which lies in the beam Bjk connecting 
i with the receptor bin j during the position k given by the 
angle 0* (see Fig. 2). Based on the geometry, the expected 
values of the Yjk are 

N 

E{Y j k) = Yl Ci 3 kXik>0 ' (1) 

For later use, let us fix some notation here. The coefficients 
Cijk give rise to two linear operators. We define C : R NS — > 
R MS by (Cx)jk = Eili CijkXik and T : R NS — > R* M5 by 
(Tx)ijk = CijkXik. We observe that, in practice, T is usually 
one to one, while C is typically not. 

III. Dynamic Tracing 

The critical part in modeling the dynamics of SPECT 
is the time behavior of activity Xi(t). If radioactive decay 
were the only significant part, a standard model of the form 
Xi(t) = Aie~ M with known decay constant A > 0 would 
apply. Dynamic SPECT is, however, based on the hypothesis 
that due to a flow in the organ, each voxel i may have its own 
individual decay profile. Therefore, more sophisticated models 
for decay with rates varying in time and space are needed. 
Based on experimental evidence for fatty acid myocardial 
viability [14], the following parametric model of decay: 

Xi(t) = ^i f c- Ait + Bie-^ 1 + d (2) 
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with unknown parameters Ai > 0, Bi > 0, C» > 0, and 
\i > 0, fXi > 0 was used in [20] and [5]. Discretizing (2) 
according to the S stops, tk = (k — l)T/S, T the total 
acquisition time, leads to a parametric model with five degrees 
of freedom for each i. 

Another approach suggested in the investigation of NMR- 
relaxation data (cf. [37]) uses a model of the form 



Xi(t) 



= I ^ 
Jo 



(A)e~ At dX. 



(3) 



Discretizing the time axis at the S stops tk now leads to a sum 
of exponentials with S degrees of freedom. 

A third approach, which we propose here, is to not insist on 
any specific form of the decay curves. That is, in each position 
i we allow for an activity profile with S degrees of freedom 
Xik, k= 1, . . . , S, assuming only that activity decays in time 



%il > %i2 > #i3 > ' * ' %iS and 



x iS > 0. 



(4). 



As compared to (2), this increases the flexibility of modeling 
decay, but increasing the degrees of freedom bares the risk 
of producing an under fitted model. Fortunately, experiments 
indicate (see Section VIII) that the model works well, and 
that its algorithmic and numeric advantages over (2) should 
be exploited. To introduce some standing notation we denote 
the set of x = {x ik } satisfying (4) by ft. 



IV. Maximum Likelihood 

Let us now derive the maximum-likelihood estimation for 
the unknown parameters x 6 ft based on the data y. Let g(y; x) 
be the probability density of the measurements Y given the 
activity curves x. Then we consider 



(ML,) 



maximize E(logg(y; x)) 
subject to x € ft 



i.e., we pick those parameters £ € ft which were the ones most 
likely to have produced the data y = {yjk}- As the Yjk are 
Poisson variables with means YliLi °ijk%ik, up to some con- 
stant the negative log-likelihood function — log#(?/;x) equals 



M s / N 
00*0 : = 53 53 ( 53 ^J* 2 ** 



■ yjk log 



i=l 



The problem of minimizing <f>(x) subject to the constraints (4) 
is feasible and admits optimal solutions since the objective 
clearly tends to infinity if at least one of the Xik — ► oo on 
the feasible set ft. On the other hand, solutions are typically 
not unique. Any solution a; € ft must satisfy the following 
Kuhn-Tucker conditions (cf. [10]). There exist multipliers 
Atfc > 0, « = 1, . . . , iV, and k = 1, ... , S such that for each 



i = 1, . . . , N we have the system (AT*) 



M 



j=l 



i'-l 



T{2 ~ 53 I CiftVfr / 53 ^'i 2 ^' 2 ] + ~ ^2 = 0 



M / 



N 



TiS-l ~ 53 I CijS-lVjS-l / 53 ^'j'S-lSi'S-l 



H-AiS-2 - XiS-i 
M / I N \ 

TiS ~ 53 I °iJSyjS I 53 °i f 3S X i'S ] + ~ ^iS 

^ik Xikfaik - ^ife+i) = 0 

k = 1, . . . , 5 - 1, XisXiS — 0 

where we used the abbreviation 

M 

nk 53^ > °' 



0 

:0 



(5) 



Clearly, any solution x has Y^i^ik^ik > 0 for fixed j,k. 
Notice also that any two solutions x,x 6 ft must satisfy 
TtiCijkXik = YsiCijkXik for every fixed j,k, since the 
functions 77 — ► 77 — yjk logr; are strictly convex for yjk > 0. As 
the operator given by (Cx)jk = J2i=i °ijk%ik is typically not 
one to one, (MLi) may have multiple solutions. Ensuring that 
C is only moderately defective is of relevance for the numerics 
and may be influenced by the discretization we choose (cf. [5]). 

Remark: We mention that the major difference with the 
widely known probabilistic model of Shepp and Vardi [34], [8] 
for static emission tomography is that, in the dynamic case, the 
model is extremely under fitted. For instance, in the situation 
of Section II, the unknown dynamic image has 64 3 degrees 
of freedom per slice, to which we fit of the order of 64 2 data 
(that is, 64 stops and 64 bins per camera cross-section). The 
situation improves if multiple camera heads are available, but 
the problem remains under fitted. To render it tractable, a prior 
model for the dynamics such as (2), (3), or (4) is required, 
represented by the constraint set ft. 

V. The EM Algorithm 

The EM algorithm is an iterative procedure to calculate 
maximum likelihood estimates. It first made its appearance 
in 1976 in a paper by Dempster et ai [12] (see also [22]) 
and its original intention was maximum-likelihood estimation 
with incomplete data. Since then, the EM algorithm has been 
applied in a much wider context, including situations in which 
the incomplete and the complete data space are defined in a 
somewhat artificial way. 

In our present situation, the maximum likelihood problem 
(MLi) involving the joint probability distribution g(y;x) of Y 
represents the incomplete data space. To introduce a complete 
data space we consider the random variables Zijk which 
present that part of the activity in voxel i that radiates toward 
the receptor bin j during stop k. These are Poisson random 
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variables with mean Cijk%ik and joint probability density 
f(z; x) depending on the parameter x G ft. 

The EM algorithm is now the following. By induction 
define a sequence of parameter estimates x a = {x? k } € ft, 
oc = 1,2, .. . according to the following rules. 

1) E step: Given x a = {xf k } € ft and the data y calculate 
the conditional expectation z a = E(Z \ Y = y^x 01 ). 

2) M step: Calculate the new estimate x a+1 € C by 
maximizing E(logf(z <x - 1 x) | y;x a ) over x e ft. 

The E step is fully explicit here. Namely, given the data y and 
the current x a we find that 



Zijk = Vjk • 



(6) 



This may, in fact, be obtained from the following. 

Lemma 1: Let Xi,.. t) X n be independent Poisson dis- 
tributed random variables with mean E(Xk) ~ Mfe- Then the 
conditional expectation 

E{(X u ...,X n )\X 1 + *-.+X n = d) =: (x u ...,x n ) 

is obtained as 

Xk = d • ~ . 

Pi + • • * + Mn 

For the M step, we maximize the log-likelihood function 
log/(2 a+1 ;x) over the x € ft. Up to constant terms, that 
entails minimizing 

A r M S 

1=1 j=l k=l 

over all x € ft. Obviously, this function is of the form 
$ — Hi^i w i tn eacn 
A/ S 

j=l k=l 

depending only on the variable x % = (xn, . . . ,#ts) 6 fli, 
the set of x x satisfying (4) for fixed i. Notice that the data 
y and the previous iterate x a do not enter directly into the 
calculation of x a+1 . Rather, is it z a which makes connection 
with the previous step. 
We have 



dxTk-^\ 

1=1 x 



Cijk = Tik 

Xik } %ik 



with Tik as in (5) and 



M 



(7) 



We now consider the maximum-likelihood problem for the 
complete data space which perforce splits into N problems 
of size 5 



(ML 2 ) 



maximize £(log/(z a ;z) | y\x a ) 
subject to x € ft. 



Setting ipik(x % ) := —Xik + Xi,k+u and with tpi defined as 
before, the ith problem becomes 



(^2,0 



minimize 



subject to ijjik{x % ) < 0 (k = 1,...,S- 1). 



Notice that in contrast with the original maximum-likelihood 
problem (MLi), we do not control the constraints xis > 0 
here, since they are built in through the objective due to 
a ts > ^* ^ n ^ act ' f° r fi xe d *> ^ consider j such that both 

> 0 and z£j k > 0. Then the corresponding term in t/\(x l ) 
works as a barrier function which forces the variable Xik to 
take strictly positive values. For the same reason, it is now 
reasonably clear that the optimal solutions x t,a+1 for (A/Z,2,t) 
are unique since the tpi are strictly convex and coercive. 

The Kuhn-Tucker conditions (cf. [10]) for the optimal 
solutions x* ,a+1 of (ML2 } %) imply the existence of multipliers 

> 0, A; = 1, . . . , S — 1, giving the system KT2,i 



Til 



Ti2 



'i2 



i\Oi _ \a 



x t,S-l 



t,S-2 



\,S-1 



= 0 



= 0 



= 0 



x. 



.ot+l 



+A* S -i = 0 



^ik X - X ifc+1 ^ifcl^^ 1 ~~ X ffc+l) — 0) 

fc = 1,2,...,S-1. 

Based on the optimality conditions (KT 2i i) t we may now enter 
a detailed analysis of the EM algorithm. 

VI. Alternating Projections 

We show that the EM algorithm presented in the previous 
section may be interpreted as an alternating projection proce- 
dure, in the sense of von Neumann, provided that the Euclidean 
metric is replaced with the Kullback-Leibler distance. The 
interested reader is referred to the Appendix for a brief account 
of von Neumann's alternating projections. 

Let n = NS and m = MS, p = NMS. The Kull- 
back-Leibler distance of a, 6 € is defined as 

p 

d(a, b) = ^2 a * lo s( a »/^) ~ a i + 

i=l 

Now d has properties which resemble those of a metric, but 
obviously lacks symmetry and is only defined for a > 0 and 
b > 0. In analogy with the orthogonal projection we define 
projection operators P£ and PJJ" associated with d. Given 
a closed convex subset A of R^. and a point x € the 
forward, respectively backward, projections of x onto A are 
defined as 



and 



PT(x) = argmincZ(a,:r) 



PT(x) = argnunti(rr, a). 



256 



IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 3, MARCH 1999 



With these notions, the M step presented in Section V turns 
out to be a Kullback-Leibler backward projection. 

Proposition 1: Let A := {v € R NMS : v ijk = Cijk^ik for 
some x e ft}. Given z<* e R NMS let i; a+1 be the backward 
projection of z a onto A, that is 

Then the solution x a+1 of the M. step satisfies = 
Ci jk x^\ for short, v a + l = Fx a + l . 

Proof: Indeed, for the proof we simply observe that with 
= Cijk%ik, for short, v = Fx, the objective function tp of 
the M-step satisfies ip(x) = J2ijk( v *jk ~ z tjk l°& v ijk) which 
up to a constant term equals d{z a ^v). Minimizing i/j(x) over 
x 6 0 is therefore equivalent to minimizing d(z a ,v) over 
v 6 A, Existence and unicity of the projection is guaranteed 
since z a > 0 and the set A contains points in R+;J* 5 . Notice 
also that A is closed, as it is the linear image of the closed 
convex and polyhedral set Q under T, and is therefore itself 
polyhedral. □ 

Based on this observation we now ask whether, similarly, 
the E step may be interpreted as a projection onto a convex set. 

Proposition 2: Let B := {z 6 R NMS : z > 0, 
J2iLi z ijk = Vjk Vj, k}. Let x a be the current EM iterate and 
let v a = IV*. Then the conditional expectation z a defined by 
the E step (6) is the Kullback-Leibler projection of v** onto 
B in the forward and the backward sense 

z" = Pb (v a ) = Pf(v a ). 

Proof: First, consider the case of the backward pro- 
jection. Clearly, z = Pjf{v a ) exists and satisfies the 
Kuhn-Tucker conditions. There exist multipliers \jk such 
that 

Zijk 

in tandem with the constraints 

N 

Summing over i for fixed j, k gives (1 + Xjk)yjk = 12i v ?jk- 
This readily implies (6). 

Next, consider the case of the forward projection z = 
Pj^(v a ). Here the Kuhn-Tucker conditions provide multipli- 
ers Xjk satisfying 

l°g(*Wv£- fc ) + Xjk = 0 V«, j, fc 

and with constraints as above. Taking exponentials and sum- 
ming over i for fixed j, k gives the same type of relation and 
again leads to (6). □ 
It remains to observe that any limit point of the sequence x a 
must be a Kuhn-Tucker point for the problem (MLx). This is a 
consequence of a general fact (see [38]), but could equally well 
be checked using the conditions (ATi) in tandem with (6) and 
(7) and the (KT 2l i). Indeed, convergence of the iterates v a ,z a 
implies convergence of x a ,a a . Formula (KKTi) implies 
k 



so the multipliers converge as well. Passing to the limit a — > 
oo therefore implies the corresponding conditions (ATi) and 
it remains to check that part of the Kuhn-Tucker conditions 
(ATi), which concerns the constraints Xi$ > 0. Clearly, 
however, if any such constraint is active Xis = 0, the 
corresponding conditions in (KTi) are satisfied. In case xis > 
0, the multiplier Xis = 0 will do. 

In consequence, we may now state the principal conver- 
gence result for the present Poisson EM algorithm. 

Theorem 1: The sequences v a = Tx a and z a generated 
by the Poisson EM algorithm based on the constraint set (4) 
converge to limit points v a — > v € A z a — ► z € B. Here, 
v is of the form v = Fx, with x € 0 a solution to the 
maximum-likelihood problem (MLi). 

Proof: First, observe that the problem (MLi) is convex 
and, due to the structure of the constraint set fi, admits 
optimal solutions. We claim that the sequences of alternating 
Kullback-Leibler projections 

converge to certain limit points v a — * v € A, z a — * z 6 B. In 
the case where A n B / 0, the alternating sequence converges 
to a common limit point v = z € AC\B. This is a known fact 
proved by Bregman (see, e.g., [2, Th. 8.1]). In case An B = 0, 
which is the more realistic one, the convergence proof is more 
complicated and only special cases have been considered. 
(See, for instance, [16], [17], [34], and [7], where the case 
of positivity constraints Xik > 0 was discussed.) Our present 
case, involving the constraint set (4), may be settled using the 
same type of reasoning. We skip over the details, which are 
tedious. 

Finally, notice that the point x satisfying v = Fx is uniquely 
determined by the injectivity of F and by closedness lies in Q. 
Inspecting the Kuhn-Tucker conditions for (MLx), we see that 
x must be a Kuhn-Tucker point and, by convexity, therefore 
solves (MLi). □ 

Remark: Consequently, the EM algorithm is best under- 
stood in the space of v variables, where it may be recast as 
an alternating Kullback-Leibler projection onto closed convex 
sets A and B, In fact, B is an afflne subspace restricted to R+, 
while A = F(Q) depends on the parameter space ^ and may 
therefore be quite complicated. The whole procedure may be 
generalized to any comparable situation with a closed convex 
parameter space Q. 

Remark: Replacing the constraint set Q by n, representing 
the nonlinear conditions (2), we obtain a nonconvex set A — 
r(fi). The EM algorithm could then still be formulated as an 
alternating projection between A and B, but with the obvious 
problems when projecting onto nonconvex sets. 

Remark: Iusem [16], among others, considers the case 
where we impose only positivity conditions Q = {x : x^ > 
0}, A = T(S7). Here, the backward projection onto A is 
particularly pleasant to calculate and leads to an explicit 
formula. 

VII. The Least Squares Approach 

In [20], the authors propose a different approach to the 
dynamic SPECT problem, which is based on the dynamic 
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model (2) and is solved using nonlinear least squares 

minimize \\Cx — y\\ 
(NLS) subject to x ik = AiC~ x ^ + Bie~^ + d 

Ai > 0, Bi > 0, d > 0, Ai > 0, 

Here, the operator C : R n — ► R m is defined as before 

N 

i=i 

Problem (NLS) was found difficult to treat in practice due to 
the highly nonlinear dependence of the X{ k on the parameters 
(cf [20], [5]). Following the ideas presented in Section III, we 
replace the parametric model (2) by (4). This leads to a linear 
least squares problem with inequality constraints 

(lls) minimize \\ Cx ~ y\\ 

^ ' subject to x6fi 

where ft is, as before, the set of x satisfying (4). Our aim now 
is to present a probabilistic model for both (NLS) and (LLS). 
This will also provide clues to possible numerical approaches 
based on a version of the EM algorithm. 

Following the pathway in Section V, we let the random 
variables Z^u and Yj k be defined accordingly, but now 
assume that the Z{j k are independent and normal variables 
with unknown mean Cijk^ik and common variance a 2 > 0. 
Consequently, the Yj k are independent normal variables with 
mean °ijkXik and variance A/a 2 . 

With these agreements, (LLS) is just the maxim urn - 
likelihood estimation problem 

, , maximize E(logg(y;x)) 
[mLl) subject to xeft 

where g{y;x) denotes the density of the normal law Y with 
mean Cx and co variance matrix Na 2 I m . Similarly, problem 
(NLS) is the maximum likelihood problem with the constraint 
set ft given by (2). 

Returning to the outline of Section V, we apply the EM 
algorithm to the present situation, letting Z, respectively Y, 
represent the complete, respectively incomplete, data spaces. 
As in the Poisson case, this generates sequences x a and 2 a , 
a = 1,2, . . . according to the following rules. 

1 ) E step. Given the current iterate x a € ft, calculate the 
conditional expectation z a = E(Z | Y — y'^x**). 

2) M step. Calculate the new iterate x a+1 € ft by maxi- 
mizing E(\ogf(z a ;x) | y;x a ) over x € ft. 

Here, }{z\x) denotes the normal law with mean Tx and 
covariance matrix a 2 I pi which is the joint density of Z. 

Proposition 3: Let A = {v € R p : v = Fx for some 
x € ft}. Then, given the result z a of the previous E step, the 
next M step reduces to calculating the orthogonal projection 

v a+l of z a omo A and taking x a+l tQ sat j sf y Px a+1 . 

Proof: Indeed, the negative log- likelihood function used 
for the M step, up to constant terms, equals \\z — Tx\\ 2 /{Na 2 ) 
and is to be minimized over x € ft. This establishes the 
statement. □ 
Let us now pass to the E step, which we would like to 
reveal as an orthogonal projection. This involves calculating 



the conditional expectation E(Z \ £ i==1 Zij k — yj^x 01 )- Due 
to independence of Z, this may be obtained from the following. 

Lemma 2: Let X = (Xi 1 . . . , X n ) be a vector of indepen- 
dcgt^tormal variables with mean E(Xi) = fa and variance 

1. Then the conditional expectation x = E(X \ Xi H 

+X n = y) satisfies 

1 1 n 

^ = -y + m 

t'=i 

More precisely, x is the orthogonal projection of . . . ,^„) 

onto the set of (xi, . . . , x n ) satisfying xi -\ h x n = y. 

Proof: In consequence of [29, Th. 2.1(viii)], a normal 
vector 

GMGM- v)) 

gives the conditional expectation 

E(b | a) ~ N(p 4- rir l (a - a), V — W+U^W). 

We apply this in the case b = (Xi,. . . , X n ) 7 a = X\-\ \-X n 

which gives the stated formula. □ 
With this observation, we find that the formula replacing (6) 
in the case of the normal variables is 

1 1 N 

Z S* = ]V^* + ^ X ?k ~nY^ ^'^ X i f k (8) 

i' = l 

and the E step is in fact realized as an orthogonal projection: 

Proposition 4: Let B = {zeR p : *ijk = Vjk Vj, A;}. 
Let x a be the current iterate generated by the previous M step 
and let v a = Tx a . Then the result z a of the next E step is the 
orthogonal projection of v a onto B. 

We are now in the position to state the convergence result for 
the EM algorithm in the case of normal variables. It is based 
on von Neumann's classical method of alternating projections. 

Theorem 2: Let A, B be defined as above and let z a , v" = 
Tx a be the sequences generated by the Gaussian EM algo- 
rithm. Then z a i v ot converge to limit points v a — * v e A, 
z a — ► z e B. Here, v is of the form v — Tx with x e ft a 
solution to the least squares problem (LLS). 

Proof: In the case where AnB ^ 0, it is well known that 
the sequence of alternating projections converges to a common 
limit point v = z 6 A n B which, by definition of A, is then 
of the form Tx. 

The more involved case occurs when AnB = 0. Following 
[3], a dichotomy appears. Either the alternating sequences 
converge to limit points u a v € A, z a — > z € B, with v, z 
realizing the distance between A,B t or there are no points 
realizing this distance, in which case the sequences tend to 
infinity \\v a \\ — ► oo, ||^ a || — ► oo, with ||v a - z a \\ approaching 
the distance between A, B. We will argue that in our situation 
the second case is impossible. 

Indeed, the distance d(A, B) not being attained implies the 
existence of a common asymptotic direction for A, B: there 
exist w ^ 0 and d € A y b 6 B having a + R+w C A, 
b + R+w C B. The second inclusion implies £^ u^* = 0 
for fixed fc, while the first gives w = Tx for some x G ft. 
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The Cijk being nonnegative, this is only possible when x = 0 
and, hence, w = 0. 

We deduce that v a — + v = Tx and z a — ► z. It is now routine 
to check that x is a Kuhn-Tucker point for the original problem 
(LLS). ^ □ 

Remark: Replacing by Q = {x : Xik > 0} gives an 
important special case. The orthogonal projection onto A = 
T(f2) is explicit and, in tandem with (8), leads to a formula for 
updating x a . We then have an iterative method for calculating 
nonnegative least squares solutions. One may compare this 
to other iterative methods for calculating nonnegative least 
squares solutions such as ISRA (see [27]). 

Let us analyze the M step a little further. Recall that up to 
a constant additive term and a constant factor, the negative 
log-likelihood function — log f(z a ,x) equals 



and is to be minimized over ft. Setting 

M M 

7ifc:=£4fc >0 and ** : = 5^**5* >0 W 



3 = 1 



3=1 



we have 

^( x ) = \ {yikxfk - 2S ik Xik) + constant terms. 

The Kuhn-Tucker conditions therefore imply the existence of 
Lagrange multipliers Xik > 0, i = 1, . . . , N and k — 1, . . . , 5 
satisfying (LKTi) 



7iS-l^tS-l ~^iS-l 
7iS^iS ~^t5 



+ Ks-2 ~ Xis-l 

+ Ais-i — A^s 



= 0 
= 0 

= 0 
= 0 



and 



Xik > ^ifc+i, A ifc (x^ - Xtfc+i) = 0, for k = 1, . . .5 - 1 

ZiS > 0, X iS XiS = 0. 

Therefore, analogously to the Poisson case, the original prob- 
lem of size NS splits into N problems of size 5. Notice that in 
contrast with the Poisson case (MLi^), we do have to control 
the constraints Xis > 0 here, which leads to the additional 
multipliers A<s. However, as we shall see, these constraints 
will generally be inactive and the multipliers will vanish. 

VIII. Numerical Approach 

In this section we shall discuss the practical aspects of both 
EM algorithms. Clearly, the algorithms may be expected to 
be slow since alternating projections are known to converge 
no better than with a linear rate. Nevertheless, the additional 
numerical stability gained may often justify the EM scheme, in 
particular, if some speed is recovered, e.g., by parallelizing the 
M step. The crucial question to be addressed before proposing 
this scheme is the following. Observe that in the dynamic 
case the M-step is intrinsically more complicated than in 



the stationary case, because it requires solving a nonlinear 
optimization problem (ML2 J i) where the static case goes with 
an explicit formula, (cf. [17], [7]). We therefore have to assure 
that a reasonably fast numerical solution is possible. In order 
to answer this question, we shall analyze the M steps of both 
the Poisson and the Gaussian model. 

The first surprising observation is that both M steps, al- 
though coming from completely different out sets, lead to 
exactly the same procedure. Recall the Kuhn-Tucker condi- 
tions (KT 2) i) for the Poisson M step presented in Section V. To 
simplify notation we omit the indexes i and a which are fixed, 
writing r k := and a k := &tk and, similarly, Xk := xf k for 
the solution and the multipliers. We then have the following. 

Lemma 3: For fixed i the problem (ML2,t) has a unique 
solution (xi,...,rrs) satisfying (4). This solution is of the 
following form: 



Xi = • - - = x ri > x n+ i 

> • ' * > #r t _!+l = ' * 



for appropriate 0 = vq < ri < • • * < r t — S with 



rCJ) 



+ • ' ' + O r 



_ 1 + 1 + ---+T r . 



(10) 



(ID 



the x^ being strictly decreasing. In particular, xi = x^ > 
^ and x s = x& < zs -. 

Proof: Uniqueness of the solution of (ML2,i) follows 
from strict convexity of the objective. Clearly, any solution, 
since it satisfies (4), has the form (10) and it only remains to 
establish (11). 

The complementarity condition in (KT2,i) implies that 
X rj = 0 for all j = 1, . . . , t, while the remaining multipliers 
may be strictly positive. Now, summing the Kuhn-Tucker 
equations (KT 2) i) in each block rj^x + 1,. .., r j separately 
gives (11). □ 

Remark: This result is interesting since it tells us that given 
any sequences a k > 0,7-* > 0, k = 1, . . . , 5, there exists a 
unique subdivision (10) such that x^\j = 1 defined 
through (11), is strictly decreasing. (Taking, for instance, 
7> = 1 we obtain the amusing observation that, given any 
sequence a k > 0, there exist a unique subdivision (10) such 
that the arithmetic means of the <Ji over each block are strictly 
decreasing and satisfy the boundary conditions x^ > <7i and 
x {t) < o-s-) 

Let us now pass to the Gaussian M step obtained in 
Section VII. Here we have the same observation which, again 
on suppressing the indexes i, a, is the following: 

Lemma 4: The solution {xu . . . , xs) of the Kuhn-Tucker 
conditions (LKTi) is unique and of the form (10), possibly 
with a different t. It admits a representation of the form 
(11) with Si replacing <j{ and 7* replacing r^. The sequence 
is strictly decreasing and, again, the limit conditions 
x\ — x^ > xs = x^ < ^- are satisfied. 

Proof: Indeed, starting out with the blocks (10), possibly 
with another £, we find again that the multipliers X rj must 
vanish. Summing the Kuhn-Tucker conditions (LKT») in each 
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broken line - model (2) 
continuous line - model (4) 



10 20 30 40 

Comparison of models (2) and (4), with 60 stops 



(a) 



broken fine - model (2) 
continuous line - model (4) 



0 20 30 40 50 

Comparison of models (2) and (4>, with 60 s 



(b) 



broken line - model (2) 
continuous line - model (4) 
dotted line - adapting (2) to (4) 



10 20 30 40 

Comparison of (2), (4) and smoothing (4) via (2) 

(C) 



broken line - model (2) 
continuous line - model (4) 
dotted line - adapting (2) to (4) 



10 20 30 40 50 60 70 

Comparison of (2). (4) and smoothing (4) via (2) 

(d) 



Fig. 3. Comparisons of models (2) and (4) using typical dynamic pixels from two simulated reconstructions on the basis of 60 (left column) and 80 (right 
column) camera positions. In the lower line an a posteriori fit of (2) to the data of (4) shows little difference with model (2). 



block would suffice, provided we knew that the final multiplier 
As, belonging to the constraint xs > 0, had vanished. Clearly, 
this is the case if xs > 0, which is consequently what remains 
to be checked. 

Assume to the contrary that the last block in (10) is zero, 
= x rt _ l +i = = xs = 0. Then the Kuhn-Tucker 
conditions give 

s 

A S = - Y, 6k 

fc=r t _i+l 

which, in view of As > 0 and 6k > 0, is only possible when 
6 rt i+1 = ... — 6s — 0. However, by construction (9) we 
have 6i > 0 so As = 0, xs > 0 and the result follows. □ 
The observation that both M steps are essentially equivalent, 
although with different data <r& > 0,t& > 0 versus 6k > 
0, 7jt > 0, gives us choice on how to perform the M step. In 
fact, the Gaussian M step uses nonnegative least squares, and 
for a moderate size S works faster than the Poisson M step 
(ML2,i). On the other hand, for a really large 5, if required, 
the Poisson M step could more conveniently be solved by an 
interior point method. In fact, a logarithmic barrier term for 



the constraints (4) leads to the objective 
s s-i 

k=i *=i 

One might perform only a few steps toward minimizing <f>^ for 
a fixed fi, > 0 before passing to another E step. Increasing the 
number of steps and reducing \i should then be controlled 
by testing decrease of the likelihood function in (MLi), 
respectively in (LLS), which, for theoretical reasons (see [12]), 
is known to decrease at each iteration if an exact EM step is 
performed. 

Experiment: Using the Poisson EM algorithm for the dy- 
namic model (4), with an M step based on a quadratic program 
(justified through Lemma 3 and Lemma 4), we reconstruct 
a typical slice of 64 x 64 pixels using 60-80 stops within 
the order of 40-50 min CPU, with up to 100 EM iterations. 
Performance can often be improved if a good prior guess is 
available. A considerable speedup could, however, be obtained 
by parallelizing the M step. A more detailed evaluation of 
our present approach including comparison to other dynamic 
methods, shall be presented in [21], 
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The experiment reported in Fig. 3 is to demonstrate that 
replacing the mixed-exponential decay model (2) by the more 
general decay model (4) usually does not lead to a loss of 
information. The four pictures show randomly chosen dynamic 
pixels from two simulated reconstructions on the basis of 60 
[Fig. 3(a) and (c)], respectively 80 [(Fig. 3(b) and (d)], camera 
positions. The plots in the top line show comparisons of typical 
decay curves produced by the models (2) (smooth semicolon 
curves) versus (4) (step-like curves). 

The Fig. 3(c) and (d) indicates that, typically, no informa- 
tion is lost on replacing (2) with (4). Namely, the smooth 
dotted lines show model (2) curves which a posteriori have 
been fitted to the step-like functions (4). Typically, these show 
little difference from the direct fitting of (2) (smooth semicolon 
type curves as before) so (2), if desired, could be retrieved 
from (4). As (4) has numerical advantages, the EM algorithm 
should in fact be built on (4). Fitting a model (2) decay curve 
to the step functions (4) may be deferred to the end of the 
procedure. 

IX. Cyclic Projections 

The fact that the Kullback-Leibler, respectively orthogonal, 
projections onto the set A cannot be calculated explicitly is 
clearly a drawback of our present approach, accounting for the 
fact that the speed is still inconvenient for clinical applications. 
In addition to parallelizing the M step, this problem may be 
addressed by the following idea. Split the projection onto A 
into two projection steps which are easier to compute and 
combine the two resulting points to approximate P^(z). More 
precisely, let us consider the sets 

Ai := {v e H p : v — Tx.xn > x i2 ,Xiz > £; 4 , . . .} 



and 



A 2 := {v € R p : v = Tx,x i2 > x^x^ > ; 



•}■ 



Then, A = A\ n A 2 and instead of projecting onto A we wish 
to use the projections onto the Ai which, as we shall see, are 
easier to perform. 

Inspecting a simple case where the sets Ai,A 2 ,B are the 
edges of an equilateral triangle shows that projecting cyclically 
onto the sets is not exactly what we want. There are, however, 
various possible ways in which the projected information onto 
the three sets Ai,A 2 ,B may be used to approach the point 
a € A, b € B, realizing the distance between the two sets. 
We propose the following scheme, which applies to both 
the Poisson and the Gaussian case with the corresponding 
interpretations. 

1) Given an iterate x a , respectively v a — rV\ do an E 
step using (6), respectively, (8). The result is z a . 

2) Replace the M step by the following. Project z a onto 
Ai, which gives v { = P^T(^ a ), i = 1,2. To form the 
new iterate x a + l , take = {v 1 + v 2 )/2, v"* 1 = 
T(x a+1 ). 

3) At the end of the procedure do a few correction steps 
by projecting v a onto the true set A, 



Let us see why this scheme, leading to explicit formulas, 
is expected to be faster than the original EM algorithm. The 
crucial observation is the following. 

Lemma 5: With the notation (7) and (5) let v 1 = -P^(^ a ) 
be the Kullback-Leibler backward projection of z a onto A x . 
Then v 1 = Fx 1 and x l is given by the following alternative. 
Either a^k-i/n^k-i < ^i^k/n^k, in which case 

1 _ i _ &i,2k-l + <7j,2fc 

x i,2k-l — x i,2k — - ~T~ 

or cri, 2 k-i/n t2 k-i > cri,2k/i~i,2k, in which case we have 

x l,2k-i = ^t,2fc-i/Ti,2fc-i) and x\ i2k = ai )2 k/n,2k- 

A similar formula is obtained for the backward projection 
v 2 = v 2 = Fx 2 , 

X. CONCLUSION 

We have presented two versions of an EM algorithm for 
dynamic SPECT where the recorded events yjk are assumed 
to be either Poisson or normally distributed. The EM algorithm 
is recognized as an alternating projection scheme in the sense 
of von Neumann for normal laws, and with respect to the 
Kullback-Leibler distance in the Poisson case. Our simulations 
indicate that a parametric model, based on experiments in 
myocardial viability studies (2) used in previous experiments 
[20], may be replaced by a less biased model (4) which, in 
addition, seems to have computational advantages. Modified 
versions of our approach (in the spirit of [28]), accounting 
for measurement noise in real clinical data, could easily be 
formulated. Related regularization techniques are discussed in 
[25] and [26]. Finally, we mention that our approach may 
be extended to a Baysian model, including prior information 
about the expected dynamic reconstructions. 

APPENDIX 

Let us briefly outline von Neumann's idea of alternating 
projections, originally formulated for hyperplanes, but later 
extended to more general convex sets. 

Suppose we are given two closed convex sets A and B 
and we are interested to find a point in the intersection 
C := AC\B, in short, a solution. Denote the orthogonal 
projections (also called nearest point mappings) onto set A, 
B by Pa, Pb and fix a starting point xq. The method of 
alternating projections generates a sequence 

Xq I— ► X\ := Pb%0 ^ %2 Pa^X > ^3 •= Pb x 2 

If C ^ 0, then the sequence (x n ) converges to some solution 
in C. Otherwise, it diverges. In the latter case, either the 
distance between the sets is attained, then the subsequences 
(#2n)> (£2n+i) converge to points realizing this distance, that 
is, x 2n — ► a € A, x 2n +i -> b € B and \\a - b\\ — dist(^, B) t 
or there are no points realizing the distance between A and 
B, in which case ||x n || — > +oc. 

The von Neumann scheme has often served as a numerical 
tool to calculate solutions to linear or nonlinear systems, in 
particular, if extended to the case of more than two sets (see 
[4] for more details and further references). 
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As the present paper highlights, the Gaussian EM algorithm 
is a particular case of von Neumann's scheme where the 
projection Pb onto B 9 the E step, has an explicit formula (8), 
while the projection Pa onto A, the M step, is more involved 
and usually requires solving a number of optimization prob- 
lems. In the static case, the M step Pa is equally convenient 
since it admits an explicit formula. Notice that in the case of 
the EM algorithm, the intersection AC\B is typically empty and 
the first part of the above dichotomy occurs, i.e., the distance 
between A and B is attained. 

In the case of Poisson distributed events, the same scheme 
is valid if the Euclidean distance is replaced with the Kull- 
back-Leibler distance (see Section IV and also [16] and [11]). 
For readers familiar with the Shepp-Vardi version of the EM 
algorithm for static emission tomography, we mention that 
in their case both projections Pjf, the E step, and PJJ~, 
the M step, are equally convenient to perform since they 
lead to explicit formulas. It is then possible to calculate the 
product PJJ" o P£~ and it is the iterative scheme based on this 
operator which, in medical imaging, is widely known as the 
EM algorithm. 
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